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The time-dependent probability distribution function of the height for the Kardar-Parisi-Zhang 
equation with sharp wedge initial conditions has been obtained recently as a convolution between 
the Gumbel distribution and a difference of two Fredholm determinants. We evaluate numerically 
this distribution over the whole time span. The crossover from the short time behavior, which is 
Gaussian, to the long time behavior, which is governed by the GUE Tracy- Widom distribution, is 
clearly visible. 



I. INTRODUCTION 

A recent experiment [l| on a film of turbulent liquid 
crystal realized the growth of a droplet made up of a sta- 
ble phase and expanding into a mctastable phase. The 
shape fluctuations of the droplet were measured with 
high precision. In approximation the top part of this 
droplet can be represented by a height function /i(x,i), 
x on the real line, and its dynamics is governed the one- 
dimensional version of the Kardar-Parisi-Zhang (KPZ) 
equation 

d t h{x,t) = ±\(d x h(x,t)) 2 + v d 2 h(x,t) + y/15 r)(x,t) (1) 

• The first term on the right hand side is a non-linearity 
which results from growth normal to the interface at lead- 
ing order in the local slope. The second term is the sur- 
face tension which tends to smoothcn the interface. The 
last term describes the random nucleation events at the 
interface and is modeled by a Gaussian white noise 77 with 
covariance 



{ri{x, t)r](x',t r )) = S(x - x')S(t - t') . 



(2) 



By scale invariance of the KPZ equation one can switch 
to dimensionlcss variables. We make the choice 



A = l, u = l/2, D = l. 



(3) 



The conversion to a general set of coefficients is discussed 
in Q, for example. We will consider the case of sharp 
wedge initial conditions, 

h(x,0) = -\x\/s with e>0,£^0. (4) 

The solution to ([T]) - (Q| is written as 

t x 2 /t\V3 



h(x, t) 



x 

24 ~ 2t 



£(x,f). 



(5) 



The flattening parabola should be viewed as the top part 
of the droplet in the experiment of [l[ and £(x, t) are the 
superimposed fluctuations. Of course, the KPZ equation 
holds in greater generality, in particular also for interface 



motion and growth models in higher dimensions. For sur- 
veys on the earlier developments we refer to Re- 
cently the KPZ equation has been used as a challenging 
test ground for nonequilibrium RG techniques [3, H| ■ 

The subtraction x 2 /2t is uniquely fixed by the require- 
ment that (£(x,i)) is independent of x for given t > 0. 
In fact, by the scale invariance of the KPZ equation, for 
fixed t also higher order correlations depend only on the 
relative distance in x. (£(x, t)) depends on t. To deter- 
mine its value is less obvious, since the construction of the 
solution to ((T|) - @ already requires a diverging uniform 
shift in the /i-direction. This is most easily explained 
for the initial condition h s (x, 0). The construction of the 
solution requires 



lim e h ' 

e->0 



(x,0) 



S(x) 



(6) 



which means that h £ (x,0) = — e _1 |x| — log(2e) with 
log(2e) diverging as e —> 0. Our convention here is to 
fix the first exponential moment of the solution to ([I} 
© as 



(e fc (*.*)) 



-x 2 /2t 



Therefore 



(exp((t/2)V3£( M ))) 



a t/24 



(7) 



(8) 



The rationale behind introducing the scale (i/2) 1 / 3 for 
£(x, t) and subtracting i/24 results from the KPZ scaling 
theory which asserts that £(x,i) = v^t + ©(t 1 / 3 ) for 
large t with Voo the true asymptotic velocity. Voo is non- 
universal, = —1/24 in our case and our units. The 
numerical factor (1/2) 1 / 3 will be discussed below. 

As in the experiment [l| , the quantity of prime physical 
interest is the fluctuating amplitude £(x,i) at fixed x, t. 
By stationarity one may pick x = and we define 



& = £(<),*). 



(9) 



To determine numerically the statistics of £t is extremely 
demanding and has never been accomplished, see Q for 
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the best results available. Thus it was a great break- 
through to have an exact formula for the distribution 
function F t (s) = Prob (£ t < s) [lMl- F t (s) can be writ- 
ten as the convolution of the Gumbel distribution with 
some function g tl 

F t (s) = l- duexp (-exp((V2) 1/3 (s-u)) g t (u) , 



(10) 

where gt turns out to be given by the difference of two 
Fredholm determinants, 

g t {u) = dot (1 - P u {B t - P Ai )P u ) - det (l - P u B t P u ) . 

The operator P u is the projection on [it, oo), while B t 
and Pai are defined in terms of the Airy function Ai 
respectively by 



dA 



Ai(z + A)Ai(z' + A) 
1 - exp(-(V2) 1 /3A) 

Ai(z - A)Ai(z' - A) 
+ l-cxp((i/2) 1 /3A) 



(12) 



and 



(z\P Ai \z') = M(z)M(i 



(13) 



We recall that the Fredholm determinant of a well- 
behaved ("trace class") integral operator with kernel A 
can be defined by 



oo 

dct(l + A) = l+£- 



(14) 



dyi ... dy m det (A{y } ;, y k )) 



j,k— l,...,m 



While (fTU)) , (fTTj) is decisive progress compared to a 
direct numerical simulation of the KPZ equation, it is 
still necessary to produce plots of the probability density 
function -F t '(s), which is the main goal of our communica- 
tion. The numerical task is demanding, since it concerns 
the computation of the Fredholm determinants, as will 
be explained in Section 2. For the data plot it is natural 
to distinguish short time, see Section 3, and long time 
behavior, see Section 4. 



II. NUMERICAL EVALUATION OF g t 

To plot the distribution F t , one needs to evaluate the 
t-dcpcndcnt family of functions gt(u), see (fTT) . Efficient, 
to some extent even optimal techniques to numerically 
compute Fredholm determinants have been developed re- 
cently [H, [HI . The basic idea is to employ a properly 
adapted discretization of the integrals appearing in the 
definition (fl4"|) of the Fredholm determinant. 

As an example, let us first consider the case of a kernel 
A(y, y') with support in [— 1, l] 2 , so that all the integrals 
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FIG. 1: (Color online) gt{u) from short times (red, larger 
amplitudes) to long times (green, lower amplitudes) for the 
values of t as listed in (1191) . 



in (fT4")l range from —1 to 1. A is assumed to be smooth 
in both variables. Then an integral is discretized as 



/i n 
- 1 t=i 



(15) 



The particular choice of the weights we and base points 
ye, £ = l,...,n, is called a quadrature rule. The sim- 
plest choice is the rectangular quadrature we = 2/rt, 
yi = — 1 + 2£/n. This choice, however, generically does 
not ensure the fast convergence of the discretized inte- 
grals. As emphasized in [15(, it is far better to use in- 
stead the Gauss-Legendre quadrature rule which leads to 
an exponentially fast convergence when n — > oo. For the 
Gauss-Legendre quadrature, ye is the £-t\\ zero of the ra- 
th Legendre polynomial P n (y) = (2™n!)- 1 9™(?/ 2 - 1)" 
and we = 2/{nP n -\(yg)P! n (yt)). The Gauss-Legendre 
quadrature is singled out as the only quadrature rule 
which gives the exact result for the integral over the func- 
tions / polynomial of degree at most In — 1. After dis- 
cretization, the Fredholm determinant reduces to a single 
n x n determinant using the von Koch formula flEl ] 



det(l + A) w 1 + V — 

Z — d rn I 



771 — 1 



ti,...,e m =i 

det(l + ^/wiA{ye 7 ye)Vwe)e,i'=i,...,n ■ 



(16) 



Since the Fredholm determinants in (|11[) involve ker- 
nels with support on R 2 , and not on [— 1,1] 2 , a change 
of variables must be performed in the integrals, using a 
function tp growing from ip(— 1) = — oo to (f(l) = oo, 



dzf(z)= / dytp'(y)f(<p(y)). 



(17) 



3 



The Frcdholm determinant is then approximated as 

det(l + A)& det (l + vWfe) (18) 
x A((p(y e ) f <p(yz>))V wi> <P' (ye)) ^, =1) ... )n . 

which converges to the true value of det(l + A) in the 
limit n — > oo. If the kernel itself is given by an integral, 
like B t of (p~2|) . this integral is also approximated by the 
Gauss-Legendre discretization (fT5j) with the same value 
of n. 

For the numerical computation of the Fredholm de- 
terminants in gt(u), we choose n = 30 and ip(y) = 
10 tan(7ry/2). All calculations are done with double pre- 
cision floating points. The values of t are chosen so to 
achieve a good resolution. In all figures we plot the 19 
values 

t = 0.25 0.35 0.5 0.75 1.2 2 3.5 6.5 13 25 50 

100 250 500 1000 2500 5000 10000 20000. (19) 

For the final integration over u in (|10[) we use an equally 
spaced grid with step size Su = 0.05 starting from u = 
and stopping after 5 consecutive values of |fft(u)| which 
are smaller than 10~ 5 . For all values of t considered, 
ranging from 0.25 to 20000, the upper bound w max is 
around u max = 3.25, while the lower bound w m in de- 
creases towards short times (u m i n = —5.65, —8, —10.9 
respectively for t = 20000, 2, 0.25). 

Since the convergence when n — > oo is slower for small 
values of u, the numerical computation of the distribution 
Ft is more difficult at short times. With n = 30, an 
accurate result cannot be achieved for times shorter than 
t « 0.25. 

The results for g t are displayed in Fig. [TJ Because 
of oscillations the function g t always takes some nega- 
tive values: it is only the convolution with the Gumbel 
distribution which yields a bona fide probability density 
function. The density F( is plotted in Figs. [2] and [4] 

III. SHORT TIME BEHAVIOR 
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FIG. 2: (Color online) Shifted and rescaled probability den- 
sity Ft (s) = aF((as + fi) from short times (red, leftmost 
curves) to long times (green, rightmost curves) for the values 
of t as listed in (|19|) . 
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FIG. 3: (Color online) First three normalized cumulants of 
the distribution Ft as a function of t. Short time expansion of 
— Ci (solid red line), C2 (solid blue line), and Cs/C^ 2 (solid 
green line) as in (|24[) and numerical evaluation of the cumu- 
lants — Ci (red circles), C2 (blue triangles), C^/C^ 2 (green 
squares) for the values of t as listed in (|19[) . The dotted lines 
represent the asymptotic Tracy- Widom values (|27|l (12911 , 



For short times a convenient choice for the fluctuating 
amplitude is 

ft(0,i) = ---logV2^+(-J Ct, (20) 

compare with (|5|) . As established in [l(| , C* converges for 
t — > to a standard Gaussian. A short time expansion 
for the cumulants of ( t is carried out in [l3jj. The short 
t numerical values are in perfect agreement with these 
results. In our display we focus on the intermediate time 
behavior for which, as before, we use the fluctuating am- 
plitude £ t with distribution F t (s). The properly scaled 
F((s) becomes Gaussian, 



lim F t (as + ju) = 




(21) 



where the parameters fj, — — 2 2 / 3 t 1 / 3 log(2-7rf) — 

2 -5/3 7r l/2 f l/6 and a 2 = 2-1/3^1/2^1/6 arc the leading 

order of the mean and variance of the distribution F t 
when t — > 0. This limit behavior is well reproduced in 
Fig.H 

From the evaluation of the distribution F t , one can 
obtain the low order cumulants. The j-th moment Mj is 
defined by 

/oo 
ds s j F{(s - 2~ 2/3 <- 1/3 log(27rf)) . (22) 
-00 

For all values of t considered, F t (s) is concentrated in the 
interval [—6,4]. Therefore numerically the moments are 
approximated by the Riemann sum of s J F((s) with step 
size Ss = 0.05 and s ranging over this interval. The three 
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FIG. 4: (Color online) Probability density function F/(s) for 
time t from short times (red, lower curves) to long times 
(green, upper curves) for the values of t as listed in (fT9l) . 
For ( -> oo, the density converges to the GUE Tracy- Widom 
distribution Fgue (upper blue curve) and, for t — > 0, F[ (s) 
becomes a Gaussian (rightmost blue curve) with mean and 
variance increasing respectively as £ -1//3 logf and t _1//6 . 



first cumulants are given by 

Ci = Mi , C 2 = M 2 -Ml , C 3 = M; :\\h\h-2M{ . 

(23) 

We compare them to the short time expansion of [13j, 
which with our conventions reads 



Ci 
C 2 

c 3 



-(t/2) 1 / 6 



(*/2) 



\3\/3 



-1/6 



3tt 

■T 

57T 



- ivs) {m 
m) m 



2/3 



(24) 



The comparison is presented in Fig. [3] The fit is fairly 
accurate at short times. The discrepancy for the smallest 
values of t signals that the discretization of the Fredholm 
determinant with n = 30 is not sufficiently accurate for 
t ps 0.25. 



IV. LONG TIME BEHAVIOR 

In the limit t — » oo, F t (s) converges to the GUE Tracy- 
Widom distribution, 



lira F t {s) = F GVE (s) = det(l - P S KP S ) , 

t— >oo 

where K is the Airy kernel 

poo 

(z\K\z') = dX Ai(z + A) Ai(z' + A) 
Jo 

_ Ai(z)Ai'(z') - Ai' (z)Ai(z') 



(25) 



(26) 



Just like in the case of a standard Gaussian, by con- 
vention the Tracy- Widom distribution is defined as in 
(|25j) . (|26|) . which thus determines the numerical factor 
(1/2) 1 / 3 . In other models the factor (1/2) 1 / 3 would have 
to be replaced by some model dependent constant. The 
normalized cumulants of .Fgue are [Itij 



lim Ci 

t— >oo 



-1.7710868074, 



lim C 2 ps 0.8131947928 , 

t— >oo 



lim C" 3/2 C 3 

t— J-oo 



0.2240842036 . 



(27) 
(28) 
(29) 



In Fig. |31 one observes how they are approached as 
t — > oo. Most significantly, the mean decays as i^ 1 / 3 , im- 
plying that the Tracy- Widom distribution is approached 
from the left. The probability density function F t '(s) is 
plotted in Fig. @] along with Fque ■ Note that the max- 
imum first swings to the left and then to the right as 
dictated by the short time behavior and also reflected by 
Cl 

The convergence to the Tracy- Widom distribution in 
the limit of long times has been established for a vari- 
ety of discrete growth models. We refer to [lTj for a 
recent survey. Fque(s) is the universal scaling function 
for droplet growth, i.e. for curved initial data. To next 
order, the difference between F t (s) and -Fgue(s) scales as 
i -1 / 3 [ill when t — > oo. The exponent —1/3 appears to 
be universal while the entire first order correction func- 
tion is model-dependent [l9| . 



V. CONCLUSIONS 

We studied numerically the one-point height distribu- 
tion function of the one-dimensional KPZ equation with 
sharp wedge initial data. Its universal part agrees very 
well with the experiment, in particular the long time con- 
vergence to the Tracy- Widom distribution and the relax- 
ation of the mean as t -1 / 3 . The scheme developed by F. 
Bornemann turns out to be highly efficient and in essence 
reduces the computation to evaluating determinants from 
a 2-parameter family of 30 x 30 matrices. In a previous 
attempt [l8T ] regular grid approximation is used which re- 
quires much larger matrices, size 80 x 80, with accurate 
results only for times longer than 10. Recently more ex- 
act solutions of the KPZ equation have become available, 
in particular the time-dependent joint height distribution 
at two spatial poi nts at the same time for sharp wedge 
initial data pit l20j. and the time-dependent height dis- 
tribution at one point for flat initial data, h(x, 0) = 
[2lj . Again the expressions involve Fredholm determi- 
nants, resp. Fredholm Pfaffians. It seems to us that the 
same numerical techniques used here could also be ap- 
plied advantageously to these other cases. 
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